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ABSTRACT 

We present radiation hydrodynamics simulations of the collapse of massive pre-stellar 
cores. We treat frequency dependent radiative feedback from stellar evolution and ac- 
cretion luminosity at a numerical resolution down to 1.27 AU. In the 2D approximation 
of axially symmetric simulations, it is possible for the first time to simulate the whole 
accretion phase (up to the end of the accretion disk epoch) for the forming massive 
star and to perform a broad scan of the parameter space. Our simulation series show 
evidently the necessity to incorporate the dust sublimation front to preserve the high 
shielding property of massive accretion disks. While confirming the upper mass limit 
of spherically symmetric accretion, our disk accretion models show a persistent high 
anisotropy of the corresponding thermal radiation field. This yields to the growth of 
the highest-mass stars ever formed in multi-dimensional radiation hydrodynamics sim- 
ulations, far beyond the upper mass limit of spherical accretion. Non-axially symmetric 
effects are not necessary to sustain accretion. The radiation pressure launches a stable 
bipolar outflow, which grows in angle with time as presumed from observations. For an 
initial mass of the pre-stellar host core of 60, 120, 240, and 480 M the masses of the 
final stars formed in our simulations add up to 28.2, 56.5, 92.6, and at least 137.2 M 
respectively. 

Subject headings: Accretion, accretion disks — Hydrodynamics — Methods: numerical 
— Radiative transfer — Stars: formation — Stars: massive 
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1. Introduction 

The understanding of massive stars still suffers from the lack of a generally accepted formation 
scenario. Despite the strong limitations in observations of massive star forming regions compared 
to their low-mass counterpart, past studies obtained common features of star formation, suggesting 
that the formation of massive stars can in first order be treated analog to low-mass star forma- 
tion. All epochs of the classical picture of star formation, such as gravitationally collapsing massive 
cores (e.g. Ho k Haschick 1986; Keto et al. 1987; Zhang k Ho 1997; Birkmann et al. 2007) and col- 
limated jets as well as wide-angle bipolar outflows (e.g. Henning et al. 2000; Zhang et al. 2002, 2007; 
Wu et al. 2005; Beuther et al. 2005a) were observed. Even circumstellar disks, e.g. IRAS 20126+4104 
(Cesaroni et al. 1997; Zhang et al. 1998; Cesaroni et al. 2005) and AFGL 490 (Harvey et al. 1979; 
Torrefies et al. 1986; Chini et al. 1991; Davis et al. 1998; Lyder et al. 1998; Schreyer et al. 2002, 
2006), or rather large scale toroid (Beltran et al. 2004, 2005) and flattened rotating structures 
(Beuther et al. 2005b; Beuther k Walsh 2008) could be revealed. Reviews of observations related 
to a proposed picture of evolutionary sequences of massive star formation are given in Beuther et al. 
(2007) or Zinnecker k Yorke (2007). 

Previous theoretical models focused on different points of view, namely the competitive accre- 
tion (Bonnell et al. 1998; Bonnell k Bate 2002; Bonnell et al. 2004; Bate 2009b,a) and the turbulent 
core model (McKee k Tan 2003), but they agree on the formation of accretion disks. 

If the formation of high-mass stars is therefore treated as a scaled-up version of low-mass 
star formation, a special feature of these high-mass proto-stars is the interaction of the accretion 
flow with the strong irradiation emitted by the newborn stars due to their short Kelvin-Helmholtz 
contraction timescale (Shu et al. 1987). Previous one-dimensional studies (e.g. Larson k Starrfield 
1971; Kahn 1974; Yorke k Kriigel 1977; Wolfire k Cassinelli 1987; Edgar k Clarke 2004) agree on 
the fact that the growing radiation pressure potentially stops and reverts the accretion flow onto a 
massive star yielding an upper mass limit of approximately 40 M . 

But this radiative impact strongly depends on the geometry of the stellar environment (Nakano 
1989). The possibility was suggested to overcome this radiation pressure barrier via the formation 
of a long-living massive circumstellar disk, which forces the generation of a strong anisotropic fea- 
ture of the thermal radiation field. Earlier investigations by Yorke k Sonnhalter (2002) tried to 
identify such an anisotropy, which they called the "flashlight effect" , in two-dimensional axially and 
midplane symmetric radiation hydrodynamics simulations, similar to our own. Their simulations 
show an early end of the disk accretion phase shortly after its formation due to strong radiation 
pressure. The final star masses are only marginally higher than the mass limit in spherical sym- 
metry, even if the frequency dependence of the radiation is considered. The radiative feedback 
was treated under the flux limited diffusion approximation (hereafter called FLD), but computed 
for several frequency bins. Due to the high computational cost of the frequency dependent FLD 
solver in Yorke k Sonnhalter (2002), it was unfortunately not possible to study the reason for the 
early fate of the disk accretion phase in detail. Krumholz et al. (2009) stated that the circumstellar 
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disk in the simulations of Yorke & Sonnhalter (2002) lost its shielding property because the disk 
region cannot be fed in axially symmetric configuration. Contrary to the stable radiation pres- 
sure driven outflows in Yorke & Sonnhalter (2002), they discovered in their own three-dimensional 
(frequency averaged) radiation hydrodynamics simulations (Krumholz et al. 2007, 2009) an insta- 
bility in the outflow region, leading to further accretion onto the disk. They propose that this 
so-called "3D radiative Rayleigh- Taylor instability" requires non-axially symmetric modes to occur 
(Krumholz et al. 2009) and that radiation pressure therefore cannot halt a flow of gas and dust 
in any direction. This explanation suffers from a lack of physical arguments and requires further 
detailed investigation due to the fact that first the classical Rayleigh- Taylor instability is a two- 
dimensional instability and secondly radiation cannot simply be treated as a fluid in general (as it 
is in the sense of a Rayleigh- Taylor instability). At least it remains unclear if this instability is the 
most important one. With the help of the unstable polar region, the most massive star in their 
simulations grows up to 41.5 Mq with an ongoing accretion phase, as the simulation is not finished 
yet. 

Contrary to this explanation of the short accretion phases in the two-dimensional simulations 
by Yorke & Sonnhalter (2002), we demonstrate here in detail the need of including the radiation 
physics at the dust sublimation front of the forming star to compute the correct anisotropy of 
the re-emitted radiation by dust grains. Due to the huge sink cells used in the simulations by 
Yorke & Sonnhalter (2002), the luminosity directly acts onto a disk region far beyond the actual 
dust sublimation front. The interaction of the radiation with the accretion flow at the dust subli- 
mation front is therefore artificially shifted to the radius of the sink cell, where the circumstellar 
disk would originally be shielded from the direct stellar irradiation. In this region the optical depth 
of the IR flux in the radial direction is much smaller than at the realistic location of the dust 
sublimation front, yielding a much higher fraction of isotropy of the radiation field. As shown in 
the well-established spherical symmetric studies, this isotropic radiation field is able to stop and 
revert the accretion flow onto the forming star. 

Our simulations focus on the accretion onto a single massive star in the center of the core. We 
incorporate the dust sublimation front of the forming star, resolve the vicinity of the star down 
to 1.27 AU, and evolve the system up to the order of 10 5 yr (ten times longer than ever studied 
before), including the whole disk accretion phase of the forming star. We scan the broad parameter 
space of numerical configurations as well as different initial conditions in several simulation series. 

In the following Sect. 2, we describe the details of the self-gravity radiation hydrodynamics code 
in use. The physical initial conditions of the pre-stellar cores as well as the numerical configuration 
of the simulations are described in Sect. 3. In Sects. 4 and 5, we present the results of one- and two- 
dimensional radiation hydrodynamics simulations of massive pre-stellar core collapses respectively, 
focusing on the radiative feedback onto the accretion flow while resolving the dust sublimation 
front. A discussion (Sect. 6) and summary of the most important results (Sect. 7) including a 
comparison of our results to the simulations by Yorke & Sonnhalter (2002) and Krumholz et al. 
(2007, 2009), a discussion of our assumptions, as well as a brief outlook on the future direction of 
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this research project, e.g. 3D simulations, close this publication. 



2. Physics and numerics 



In this section, we outline the ingredients and default numerical configuration of the self-gravity 
radiation hydrodynamics code we use to model the collapse of massive pre-stellar cores. The first 
Sect. 2.1 comprises the motivation for our choice of a grid in spherical coordinates and highlights 
the step forward in resolution we obtain in our simulations compared to previous research. The 
following two Sects. 2.2 and 2.3 describe the features and the configuration of the hydrodynamics 
solver including full tensor viscosity, for which we use the open source magneto-hydrodynamics 
code Pluto3 (Mignone et al. 2007). Further sections describe our newly developed modules for 
self-gravity (Sect. 2.4) and frequency dependent approximate radiation transport (Sect. 2.5). We 
close this section with the description of the pre-calculated, tabulated dust (Sect. 2.6) and stellar 
evolution model (Sect. 2.7) used in the simulations. 



In our simulations, we use a time independent grid in spherical coordinates with logarithmically 
increasing radial resolution towards the center. The usage of a radially increasing resolution towards 
the center guarantees the possibility to study the radiative feedback in the central core regions 
down to a minimum grid cell size of Ar x r A9 =1.27 AU x 1.04 AU. An example of such a two- 
dimensional grid is displayed in Fig. 1. This type of grid is well adapted for the analysis of the 
interaction of an accretion flow onto a massive star along with the stellar irradiation it generates, 
because the stellar gravity as well as the stellar radiative force are aligned with the radial coordinate 
axis. Furthermore, in contrast to e.g. cartesian coordinates, the usage of spherical coordinates 
guarantees a strict angular momentum conservation. The polar discretization A9 of the grid is 
uniformly fixed and covers an angle of tt/2 from the top polar axis to the forming disk midplane, 
assuming midplane symmetry as in Yorke & Sonnhalter (2002). The polar resolution r A9 of the 
spherical grid automatically increases towards the regions of interest around the centrally forming 
star, where high resolution is desired. To even enhance this focus on the inner parts of the pre- 
stellar core and saving computational effort in the outer parts far away from the dust radiation 
interaction layer, we choose a logarithmically increasing radial resolution of the grid. Thus the 
non-adaptive grid setup impeded the study of potential fragmentation in the outer core regions. 
The radial resolution Ar at a radius r of a computational domain with N r grid cells in the radial 
direction is given by 



with / = log (r max /r min ) /N r , where 

r mm and r max represent the inner and outer radius of the 
computational domain. A comparison of our achieved resolution to previous massive pre-stellar 



2.1. Discretization of the computational domain 
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(a) Global image of the total computational domain up to (b) Zoom-in image of the central 100 AU x 100 AU. 
the outer radius of r max = 0.1 pc. The innermost cells have a resolution of 

Ar x r A6 =1.27 AU x 1.04 AU. 

Fig. 1. — : Two-dimensional grid (64 x 16) in spherical coordinates with logarithmically increasing 
radial resolution, a central sink cell of radius r m ; n = 10 AU, and an outer boundary at r max = 0.1 pc. 
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core collapse simulations by Yorke &: Sonnhalter (2002) and Krumholz et al. (2007, 2009) is given 
in Table 1. 

The forming high-mass proto-stellar object at the center of the core is represented by a ded- 
icated stellar evolution model (presented in Sect. 2.7) inside the central sink cell with radius r m \ n 
at the origin of the coordinate system using pre-calculated stellar evolutionary tracks for accreting 
high-mass proto-stars (Hosokawa &: Omukai 2009). 



2.2. Hydrodynamics 

To follow the motion of the gas, we solve the equations of compressible hydrodynamics 

3t p + V • (p u) = (2) 
d t (pu) + V(puu + P) = pa (3) 
dt E + V- ((E + P) u) = pu-a-VF tot (4) 

with the acceleration source term a = X^i^i which includes the additionally considered physics 
to the equations of gas dynamics (Euler equations) such as shear viscosity (ai), central gravity 
of the forming star (0*2), self-gravity (03), and radiation transport and stellar radiative feedback 
(04). -Ftot denotes the flux of the total radiation energy density. These additional components are 
described in the following subsections. The evolution of the gas density p, velocity u, pressure P, 
and total energy density E is computed using the open source magneto-hydrodynamics code Pluto3 
(Mignone et al. 2007). 

Pluto is a high-order Godunov solver code, i.e. it uses a shock capturing Riemann solver within 
a conservative finite volume scheme. The numerical configuration of our simulations makes use of a 
Strang operator splitting scheme for the different dimensions (Strang 1968). Our default configura- 
tion consists further of a Harten-Lax-Van Leer (hll) Riemann solver and a so-called 'minmod' flux 
limiter using piecewise linear interpolation (plm) and a Runge-Kutta 2 (RK2) time integration, also 
known as the predictor-corrector- method, compare van Leer (1979). Therefore the total difference 
scheme is accurate to 2 nd order in time and space. 

To close the system of Eqs. (2) to (4), we use an ideal gas equation of state 

P = ( 7 -l)£ int , (5) 

which relates the gas pressure P to the internal energy £j nt = E — O.bpu 2 . The adiabatic index 7 
is set to 5/3. 

To limit the range of densities, the so-called floor value of the density is chosen to be po = 
10~ 21 g cm -3 . This floor value occurs during the simulations only in regions where the radiation 
pressure driven outflow is depleting the density of the corresponding grid cells in the radially 
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Table 1. Resolution of different radiation hydrodynamics simulations of a collapse of a slowly 

rotating massive pre-stellar core 



Authors 


Resolution in 
lowest resolution 


AU in regions of 
highest resolution 


Radius of 
sink cells (AU) 


Yorke & Sonnhalter (2002) 


320 2 


80 2 


80 


Krumholz et al. (2007) 


966 3 


7.5 3 


- 30 


Krumholz et al. (2009) 


645 3 


10 3 


- 40 


This study ID 


1540 


0.08 


1.0 


This study 2D 


2319 x 1911 


1.27 x 1.04 


10.0 



Note. — The simulations of Yorke &; Sonnhalter (2002) were performed on a non-adaptive two- 
dimensional grid in cylindrical coordinates with three levels of refinement. The given resolution 
(Ar x Az) of Yorke & Sonnhalter (2002) represents the case of a M core = 60 M pre-stellar core. 
The resolution for the lower mass M corc = 30 M Q collapse was a factor of two better. The reso- 
lution for the higher mass M core = 120 M collapse was a factor of two worse. The simulations 
of Krumholz et al. (2007, 2009) were performed on a three-dimensional cartesian adaptive mesh 
refinement (AMR) grid. The given resolution (Ax x Ay x Az) represents the lowest and highest 
resolution during the simulation. The radiation from the sink cells is added to their computational 
grid using a smooth weighting function inside the so-called accretion radius, which is four times 
the highest resolution of the grid. (Krumholz 2005). The resolution of our own grids in spherical 
coordinates is given in units of arc length (Ar x (r AO)). 
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outward direction. Thus, the choice of the floor value does not influence the level of accretion onto 
the newly forming star we investigate. 

The various sources of additional acceleration a = £^ a, that enter the equations of hydrody- 
namics in Eqs. (3) and (4) are discussed in the following sections and include viscosity, gravity of 
the central star as well as self-gravity of the core, and radiative feedback. 



2.3. Viscosity 

In the two-dimensional simulations we consider physical shear viscosity of the circumstellar 
disk medium to mimic the effect of angular momentum transport (via e.g. the magneto-rotational 
instability, spiral arms, disk winds and jets). Two-dimensional axially symmetric simulations with- 
out any shear viscosity yield the formation of ring instabilities in the circumstellar disk. The rings 
would be unstable if non-axially symmetric modes were allowed, leading to the formation of spiral 
arms and therefore to angular momentum transport as discussed by Yorke et al. (1995). 

Full physical tensor viscosity is included in the current version of the open source magneto- 
hydrodynamics code Pluto3 (Mignone et al. 2007). Viscosity enters the equations of hydrodynamics 
in Eqs. (3) and (4) as an additional source of acceleration 

oi = vn. (6) 

The components of the viscous stress tensor II are given (in cartesian coordinates) by 

Rij = V (^djUi + diUj - ^8ijd k uk \ + Vb^ijOkUk (7) 

with the shear viscosity 77, the bulk viscosity 77b, and the Kronecker symbol Sij. Further details on 
the analytical treatment of viscosity can e.g. be found in Landau & Lifshitz (1987). We assume for 
the bulk viscosity 77b = 0. 
The (shear) viscosity 

rj = pv (8) 

is described via the so-called a-parameterization of Shakura &: Sunyaev (1973), in which the dy- 
namical viscosity v is set proportional to the product of a typical velocity and length scale of the 
system under investigation, here the local sound speed c s and pressure scale height H: 

v = a c s H (9) 

We further approximate the local pressure scale height H by 



H = ^ (10) 



with the keplerian angular velocity 



<*r)-J™«l (11) 
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derived from the equilibrium between gravity and the centrifugal force. The mass M(r) inside the 
radius r is calculated by the spatial integral of the density distribution plus the central stellar mass 
M* inside the sink cell: 

i*r pit 

M(r) = M* + 2vr dr d9 p(r,0) r 2 sin(0). (12) 
Jo Jo 

Using the relation (10) we substitute the local sound speed in Eq. (9) yielding 

v = a n K (r) H 2 . (13) 
Introducing the dimensionless parameter H/R, the aspect ratio of the circumstellar disk, leads to 

v = a n K (r) R 2 (^) (14) 

with the cylindrical radius R = r sin(#). If the viscosity is an effect of turbulent trans- 
port of angular momentum, e.g. by the magneto-rotational instability Balbus & Hawley (1991); 
Hawley & Balbus (1991); Balbus (2003) or the baroclynic instability Klahr & Bodenheimer (2003), 
it is observed that the strength of the stresses is proportional to the thermal pressure. This is the 
fundamental assumption of the a-Ansatz by Shakura & Sunyaev (1973). This relation (Eq. (9)) 
holds because hotter and thicker disks can support higher levels of turbulence. The situation is 
reversed for self-gravitating disks. Here, hot disks are usually Toomre stable and will not produce 
gravito-turbulence. On the contrary, the disks will cool down to the marginally unstable Toomre 
values and establish a turbulent state where the level of turbulence is set by the equilibrium of 
energy release and radiative cooling (Gammie 2001). For that reason, we choose a viscosity pre- 
scription independent on the actual disk temperature (e.g. a fixed H/R ratio of 0.1 and a = 0.3) 
but only on the local mean (keplerian) rotation profile. This way, we ensure that cool and thin 
disks can obtain the high viscosity values they deserve. Our Ansatz is equivalent to the so-called 
/3- viscosity Ansatz for self-gravitating disks by Duschl et al. (2000), which is also independent on 
temperature, with the /3-parameter of /3 ~ 3 * 10~ 3 . 



2.4. Gravity 

The calculation of the gravitational potential $ is split into the gravity of the central star in 
the sink cell <3?* and the self-gravity of the mass in the computational domain <I> S g: 

$ = + $ sg (15) 

The associated accelerations c?2 + 03 enter the hydrodynamics as a source term for momentum and 
energy in Eqs. (3) and (4). 

The acceleration vector 02 of the gravity of the central star is given analytically by 

- , -CM* G M*_, GM^ 
CL2 = -V <P* = V = o r e r = ^—e r . (16) 

£ 
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Such external gravity (from point sources) is supported in Pluto3 by denning the gravitational 
potential or the resulting acceleration vector 3,2- 
The acceleration 0,3 due to self-gravity is given by 

a 3 = -V $ sg , (17) 

in which the gravitational potential & sg is determined via Poisson's equation: 

V 2 $ sg = 4tt G p. (18) 

We implemented a solver of Poisson's equation into our version of the Pluto code in a modular 
fashion. The module solves the equation via a diffusion ansatz. The desired approximate matrix 
inversion is done using the so-called GMRES method. The accuracy of the Poisson solver, i.e. the 
abort criterion for the approximate matrix inversion, is choosen to 0.001% relative accuracy of the 
gravitational potential (A<I> S g/<]? S g < 10~ 5 ). 

The outer radial boundary values of the gravitational potential are calculated via a Taylor 
expansion of the density distribution, described for example in Black & Bodenheimer (1975). Sev- 
eral tests we performed indicate that it is sufficient to just account for the monopole solution of 
the Taylor expansion, i.e. the total mass of the core. In our default configuration of the pre-stellar 
cores (see Sect. 3), the mass distribution is perfectly spherically symmetric at the beginning of the 
simulation and afterwards becomes highly concentrated in the inner region of the computational do- 
main far away from the outer boundary, both yield analytically the monopole solution at the outer 
boundary. To control the resolution, which is necessary to resolve the physics of self-gravity cor- 
rectly, e.g. preventing artificial fragmentation, we monitor the so-called Truelove criterion, derived 
in Truelove et al. (1997). The criterion requires to resolve the Jeans length 




at least by a priorily defined number of grid cells. The inverse of the number of necessary grid cells 
per Jeans length is the so-called Jeans number 

J = (20) 

Truelove et al. (1997) suggested at least a Jeans number of J < 0.25. The worst case during 
our simulations occurs in the high-density region around the forming star with approximately a 
maximum density of p < 10~ n g cm -3 , a minimum temperature of T ~ 100 K and a resolution of 
the order of Ar « 1 AU. This leads to a Jeans number of J « 0.09, i.e. the Jeans length Aj is at 
least resolved by 11 grid cells. 



2.5. Radiation transport 



The importance of the frequency dependence of the stellar spectrum when calculating the 
radiative feedback of a massive star was already shown in radiation hydrodynamics studies by 
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Yorke & Sonnhalter (2002) and Edgar & Clarke (2003) as well as in radiation transport simulations 
without hydrodynamics by Krumholz et al. (2005). On the other hand, no frequency dependent 
radiation hydrodynamics study related to massive star formation was carried out since the work by 
Yorke &c Sonnhalter (2002) in more than one dimension and due to the huge computational overhead 
of their frequency dependent FLD routine it was neither possible to study a large number of different 
initial conditions (to scan the parameter space), nor to perform high-resolution simulations of the 
accretion process. 

Furthermore, stellar radiative feedback onto the dynamics of the environment plays a crucial 
role in the formation of massive stars. On the one hand, the heating will probably prevent further 
fragmentation of the cloud by enhancing the Jeans mass (e.g. Krumholz et al. 2007). On the 
other hand, the dusty environment feels the radiative force when absorbing the radiation due to 
momentum conservation, which potentially stops the accretion process for highly luminous massive 
stars. 

To study the radiative feedback of massive stars on their own accretion stream in one-, two-, 
and three-dimensional simulations we implemented a fast, robust, and accurate frequency depen- 
dent radiation transport solver in spherical coordinates into our version of the Pluto code. To 
achieve a fast solver for the frequency dependent problem we split the radiation field into the stel- 
lar irradiation and thermal dust emission. The basic methodology of the hybrid scheme is to 
perform a frequency dependent ray-tracing step for the stellar irradiation and shift the re-emission 
of the photons by the dusty environment to a frequency averaged flux limited diffusion solver in 
the equilibrium temperature approximation. A derivation of the hybrid scheme and numerical 
details of the implementation of this newly developed radiation transport method are given in 
Kuiper et al. (2010), including a detailed comparison of the method with the standard radiation 
transport benchmark test by Pascucci et al. (2004). 



2.6. Dust model 

For the implementation of realistic mass absorption coefficients k {v) for the frequency depen- 
dent radiation transport module, we use an opacity table of Laor k, Draine (1993), including 79 
frequency bins, shown in Fig. 2. The opacity table covers the full frequency range from microwave 
and infrared radiation up to soft x-rays. It describes a mixture of dust grains in the size range 
between 0.005 to 10.0 (im.. The grains are taken to be spherical and consist out of amorphous sili- 
cate with a composition like that of olivine. The model does not include ice mantles, which will 
roughly double the opacity at temperatures below ~ 100 K (e.g. Ossenkopf & Henning 1994). As 
shown in Fig. 2 this dust grain mixture takes into account the strong absorption/emission features 
at 9.7 |im and 18 fim observed in the interstellar medium. The corresponding frequency averaged 
Planck- and Rosseland mean opacities are shown in Fig. 3 as a function of temperature. 

Aside from the dust opacities, the opacity of a given grid cell depends also linearly on the local 
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Fig. 2. — : Frequency dependent mass absorption coefficients in tabulated form from 

Laor k Draine (1993). 
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Fig. 3. — : Calculated Rosseland and Planck up mean opacities as a function of dust tempera- 
ture. Here, the possible evaporation of dust grains at high temperatures (and/or low densities) is 
neglected, but considered in the dust to gas mass ratio Md us t/M gas of each grid cell, cp. Fig. 4. 



dust to gas mass ratio. The initial dust to gas mass ratio (Afd us t/M gas ) is fixed to 1%. Gas and dust 
is treated as a single fluid, so the dust to gas mass ratio only shrinks due to possible evaporation of 
the dust grains in hot regions (around the central massive star). The local evaporation temperature 
of the dust grains is calculated by using the formula of Isella & Natta (2005) 

T cvap = g P p (21) 

with g = 2000 K, f3 = 0.0195, and the gas density p given in g cm~ 3 . The formula describes 
a power-law approximation to the evaporation temperatures T evap determined by Pollack et al. 
(1994). A smooth spatial and temporal transition of the associated dust to gas mass ratio between 
completely evaporated and condensated regions is achieved via the transition function 



M gas V ^gas J V ' V 100 

The transition slope is displayed in Fig. 4 as a function of the temperature for a high gas density 
of p = 10~ 10 g cm" 3 as well as for the floor value of the density po = 10 -21 g cm" 3 . 



-14- 



1.0 




0.0 1 L ■ ■ ■ 7—7-TT === — i 

500 1000 1500 2000 2500 3000 
T[K] 

Fig. 4. — : Transition slope of the local dust to gas mass ratio as a function of the temperature due to 
evaporation of dust grains for two different gas densities. The vertical lines mark the corresponding 
evaporation temperatures. 



2.7. Stellar evolution model 

The evolution of the central star, described by a central sink cell, is coupled to the hydrody- 
namics of the pre-stellar core by measuring the mass flux into the sink cell. The initial mass of 
the star at the beginning of the simulation is simply given by the integral over the initial density 
distribution up to the radius r m j n of the sink cell and is therefore in all cases less than a few percent 
of 1 M.q. The mass, which enters the sink cell during the hydrodynamics is assumed to be accreted 
onto the central star. From the mass flux pu into the sink cell during the timestep At we calculate 
the accretion rate M onto the central star via 

M = 2vr / d6 pu-e r r^ in sin 6. (23) 
J o 

Integrating the accretion rate M over the time yields the growth of the stellar mass M* : 

/t+At 
Mdt = M*(t) + MAt. (24) 

A potential decrease of the accretion rate due to outflows or jets (driven by magnetic forces) in 
the inner region below the sink cell radius, is currently not considered. The total luminosity L to t 
is given by the sum of the accretion luminosity L acc and the stellar luminosity L*: 



£tot(*) = L acc(t) + £*(*)• 



(25) 
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The accretion luminosity is directly calculated from the hydrodynamics simulation via 

^acc = (26) 

with the stellar radius i?* . The stellar luminosity and the stellar radius are obtained via fits to the 
pre-calculated evolutionary tracks by Hosokawa & Omukai (2009). These evolutionary tracks of 
massive stars depend on the stellar mass as well as on the actual accretion rate. We use polynomial 
fits to the mass relation up to 10 th order for separated mass ranges (an example of these fits is 
shown in Fig. 5) and linear regression for the dependency on the instantaneous accretion rate. 
For stellar masses below the accessible data (0.05 M Q in the worst case) the stellar luminosity is 
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Fig. 5. — : Polynomial fits to the stellar luminosity as a function of the stellar mass as calculated by 
Hosokawa & Omukai (2009). The data points represent an evolving massive star with an accretion 
rate of 10 -3 Mq yr _1 . The mass range was split into two regimes above and below 5.5 M Q (at the 
sharp bend) and each part is fitted by a polynomial up to 10 th order (solid lines). 



assumed to be negligible and the stellar radius is assumed to be constant up to the first data point. 
Given the stellar radius and total luminosity, the stellar effective temperature T* is calculated 

from 

Ltat = 4vr <j S b T*. (27) 
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3. Initial conditions and numerical configuration 

Using the newly developed modules of the self-gravity radiation hydrodynamics code presented 
in Sect. 2, we performed multiple simulations of collapsing massive pre-stellar cores. Most of the 
simulations were performed either to scan the huge numerical parameter space of the setup to 
guarantee significant results or to explore individual physical initial conditions. An overview of 
the 22 simulations evaluated is presented in Table 2 and Table 3 for one- and two-dimensional 
simulations respectively. 

Aside from varying one specific parameter of the initial condition or the numerical configuration 
in each simulation series, most of the initial conditions and the physics considered in the simulations 
stay the same. 

Our basic initial condition is very similar to the one used by Yorke &; Sonnhalter (2002). We 
start from a cold (To = 20 K) pre-stellar core of gas and dust. The initial dust to gas mass ratio 
is chosen to be Md ust /M gas = 1%. The model describes a so-called quiescent collapse scenario 
without turbulent motion (u r = uq = 0). In non-spherically symmetric two-dimensional runs the 
core is initially in slow rigid rotation (ju^/R = Qq = 5 * 10~ 13 Hz). The rotation speed of f^o 
results roughly in an equilibrium between gravity and centrifugal force at the outer core radius 
r max in the case of the lowest mass core of M core = 60 M . The outer radius of the cores is fixed 
to r max = 0.1 pc and the total mass M core varies in the simulations from 60 up to 480 Mq. The 
initial density slope drops with r~ 2 . The highest mass case of 480 M with a mean density of 
p ~ 8 * 10~ 18 g cm" 3 («2* 10 6 cm -3 ) denotes the upper mass limit of such a pre-stellar core we 
would expect from observations. A brief overview of these physical initial conditions of the massive 
pre-stellar cores studied is given in Table 4. 

The simulations are performed on a time independent grid in spherical coordinates (see Sect. 2.1). 
The radially inner boundary of the computational domain is a semi-permeable wall towards the 
forming star, i.e. the gas can enter the central sink cell, but it cannot leave. The outer radial 
boundary is a semi-permeable wall as well. The mass can be pushed out of the computational 
domain (by radiative forces) but no mass is allowed to enter the computational domain. The semi- 
permeable outer boundary implies the assumption that the collapsing core is mostly isolated from 
its large-scale environment. This limits the extent of the potential mass reservoir for the forming 
massive star to the initially fixed mass of the pre-stellar core M core . 

The remaining numerical parameters are determined in several simulation series. The reso- 
lution of the computational domain, which is necessary to follow the radiation and fluid physics 
as well as its interactions, is determined in several so-called convergence runs, see Sects. A.l and 
A. 2 for non-rotating and rotating cores respectively. The highest resolution of the non-uniform 
grid is chosen around the forming massive star, afterwards the resolution decreases logarithmically 
in the radial outward direction. The default resolution goes down to (Ar) m i n = 0.08 AU and 
(Ar x r A#) min = 1.27 AU x 1.04 AU for the one- and two-dimensional simulations respectively. 
The accurate size of the sink cell is determined in parameter scans presented in Sects. 4.1 and 5.1 for 
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Table 2. Overview of spherically symmetric massive pre-stellar core collapse simulations 









Resol. 


^min 


•^core 


ts 


^end 


Label 


Dim. 


Grid cells 


(AU) 


(AU) 


(M©) 


(kyr) 


(kyr) 


ID Convergence runs 




Sect. A.l 












lD-Convergence32 


ID 


32 


0.36 


1 


60 


67.6 


310** 


1 D- Conver gence64 


ID 


64 


0.17 


1 


60 


67.6 


163** 


lD-Convergencel28 


see 'lD-Mcore60Msol' 










lD-Convergence256 


ID 


256 


0.04 


1 


60 


67.6 


188* 


ID rmin runs 




Sect. 4.1 












lD-rminlAU 


ID 


99 + 128 


1.0 


1 


60 


67.6 


204* 


lD-rmin5AU 


ID 


95 + 128 


1.0 


5 


60 


67.6 


282** 


lD-rminlOAU 


ID 


90 + 128 


1.0 


10 


60 


67.6 


283** 


lD-rmin80AU 


ID 


20 + 128 


1.0 


80 


60 


67.6 


293* 


ID Mcore runs 




Sect. 4.2 












lD-Mcore60Msol 


ID 


128 


0.08 


1 


60 


67.6 


218** 


lD-Mcorel20Msol 


ID 


128 


0.08 


1 


120 


47.8 


54** 


lD-Mcore240Msol 


ID 


128 


0.08 


1 


240 


33.8 


39** 


lD-Mcore480Msol 


ID 


128 


0.08 


1 


480 


23.9 


10** 



Note. The table is structured in blocks of topics and their corresponding sections. For 
each run the label, the dimension, the number of grid cells, the resolution of the best resolved 
region around the central star, the radius r m i n of the central sink cell, the initial mass M core of the 
pre-stellar core, its corresponding free fall time tg = tt r^lax/V^ G M core , and the period t cnc j of 
evolution simulated are given. A '*' in the i cn d column denotes that the whole accretion phase of 
the star has been computed. A '**' denotes that the computation has been stopped at the point 
in time when no mass is left in the computational domain. 
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Table 3. Overview of axially symmetric massive pre-stellar core collapse simulations 









Resol. 


^min 


-'"core 


iff 


^end 


Label 


Dim. 


Grid cells 


(AU) 


(AU) 


(M ) 


(kyr) 


(kyr) 


2D Convergence runs 




Sect. A.2 












2D-Convergence32xl6 


2D 


32 x 16 


2.69 x 1.11 


10 


60 


67.6 


62 


2D-Convergence64x4 


2D 


64 x 4 


1.27 x 4.18 


10 


60 


67.6 


93 


2D-Convergence64x8 


2D 


64 x 8 


1.27 x 2.09 


10 


60 


67.6 


691** 


2D-Convergence64xl6 


see '2D-Mcore60MsoP 










2D-Convergencel28x32 


2D 


128 x 32 


0.61 x 0.51 


10 


60 


67.6 


33+ 


2D rmin runs 




Sect. 5.1 












2D-rmin80AU 


2D 


64 x 16 


7.25 x 8.21 


80 


60 


67.6 


251* 


2D-rminlOAU 


see '2D-Mcore60Msol' 










2D-rmin5AU 


2D 


64 x 16 


0.69 x 0.52 


5 


60 


67.6 


631+ 


2D-rminlAU 


2D 


64 x 16 


0.17 x 0.11 


1 


60 


67.6 


92+ 


2D Mcore runs 




Sect. 5.2 












2D-Mcore60Msol 


2D 


64 x 16 


1.27 x 1.04 


10 


60 


67.6 


939** 


2D-Mcorel20Msol 


2D 


64 x 16 


1.27 x 1.04 


10 


120 


47.8 


489** 


2D-Mcore240Msol 


2D 


64 x 16 


1.27 x 1.04 


10 


240 


33.8 


226** 


2D-Mcore480Msol 


2D 


64 x 16 


1.27 x 1.04 


10 


480 


23.9 


41+ 



Note. The table is structured in the same way as Table 2. Simulations, which are still 
running, are marked by an additional '+' in the £ e nd column. 
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non-rotating and rotating cores respectively. While we started the first two-dimensional simulations 
with a radius of the inner sink cell of r m j n = 80 AU analog to the case of 'F60' in Yorke & Sonnhalter 
(2002) we experienced that it is necessary to shrink the size of the sink cell down to a value smaller 
than the distance from the dust sublimation front to the forming massive star, at least from the 
point in time on at which the radiative force becomes a serious counterpart to the gravity. Oth- 
erwise a huge sink cell give rise to an artificially high radiative feedback and therefore limits the 
stellar mass reached in the simulations dramatically. The default radius of the sink cell is chosen 
to be r m i n = 1 AU and r m i n = 10 AU for the one- and two-dimensional simulations respectively. In 
axially symmetric (two-dimensional) runs, physical shear viscosity is used to maintain the accretion 
flow through the growing circumstellar disk. Therefore, we adopted the well-known a-parametri- 
zation model for shear viscosity of standard disk theory (Shakura & Sunyaev 1973). We performed 
several simulations with varying normalization values for the physical a-viscosity, which yield the 
formation of a stable accretion disk for a range of a-values from 0.1 up to 1.0. Apart from these 
runs the normalization of the viscosity was fixed to be a = 0.3 here. A theoretical estimation on 
the a-values of massive accretion disks was presented in Vaidya et al. (2009). 

In previous test runs, we studied several non-radiative and radiative physics. We performed 
isothermal and adiabatic collapse simulations as well as gray and frequency dependent radiation 
transport with and without radiation pressure feedback from the star or the diffuse thermal radi- 
ation field. In this paper, we confine ourselves to present only the most realistic runs including 
frequency dependent radiation transport as well as full radiative feedback. 

4. Spherically symmetric accretion 

As mentioned in the introduction, spherically symmetric accretion onto a massive star is po- 
tentially stopped by its growing radiation pressure. First, we present our results of one-dimensional 
simulations used to fix the numerical parameters of the setup, namely, the resolution of the compu- 
tational grid (see Sect. A.l) and the radius of the inner sink cell (Sect. 4.1). Afterwards, we analyze 
simulations with varying initial core masses M core to determine the upper mass limit, if such a limit 
exists, for our specific model (the chosen dust and stellar evolution model, the configuration of the 
hydrodynamics solver as well as the treatment of radiation transport). To re-perform these one- 
dimensional simulations on our own instead of just referring to Larson Sz Starrfield (1971), Kahn 
(1974), Yorke & Kriigel (1977), and Wolfire & Cassinelli (1987) allows us to directly compare the 
results found for spherically symmetric accretion flows with subsequent simulation results in higher 
dimensions. 
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4.1. Parameter scan of the size of the sink cell: The influence of the dust 

sublimation front 

4-1.1. Simulations 

In order to limit the run time of the simulations to an adequate amount, the formation and 
evolution of the central proto-star cannot be included in the computational domain. In fact, the 
radially inner computational boundary defines the radius of a so-called sink cell. The mass flux 
into this sink cell defines the accretion rate onto the proto-star, which is assumed to form in the 
center of the pre-stellar core. Inside of this sink cell the stellar properties such as luminosity and 
radius are taken from pre-calculated stellar evolutionary tracks. We use therefore recent results 
for the evolution of accreting high- mass stars by Hosokawa & Omukai (2009). In the following, we 
study the influence of the location r m \ Q of this inner boundary on the resulting accretion rate onto 
the evolving massive star. We check this dependency in four simulations with a radius of the inner 
sink cell of r m ; n = 1, 5, 10, and 80 AU. To decouple the results from the dependence on resolution 
(see Sect. A.l) the size of the grid cells was fixed to be Ar = 1 AU up to a radius of 100 AU. So the 
different simulations use 99, 95, 90, and 20 grid cells up to 100 AU respectively. Behind this inner 
region, the grid resolution decreases logarithmically throughout additional 128 grid cells from 100 
AU up to 0.1 pc. The initial conditions and numerical parameters of these runs are described in 
Sect. 3 and the simulations are performed for an initial core mass of M core = 60 M Q . We follow the 
long-term evolution of the runs for at least 200 kyr, representing 3.0 free fall times. The resulting 
accretion flow onto the forming star as well as the deviations of the simulations from the run with 
the smallest sink cell r m i n = 1 AU are displayed in Fig. 6. 

4-1-2. Conclusions 

The first absorption of stellar irradiation takes place directly behind the dust sublimation 
radius r su bi- If the radius of the central sink cell r m j n exceeds this dust sublimation radius, this 
interaction is artificially shifted to r m ; n . Due to the fact that the generalized Eddington limit is 
independent of the radius (the stellar gravity and the stellar radiative flux both drop with r -2 ), 
the shift of this first transfer of momentum from the stellar irradiation to the dust flow should be 
independent of the radius of the sink cell. Secondly, the absorption of stellar irradiation heats up 
the region behind the dust sublimation radius respectively. The thermal radiative flux from this 
region outwards slows down the gravitationally in-falling accretion flow. In general, this interaction 
depends on the radius, which defines the temperature of the heated region, the velocity of the 
accretion flow and the opacity of the corresponding dust. 

In the plot of the resulting accretion rates (Fig. 6, upper panel) only slight deviations of the 
run with the largest sink cell radius r m ; n = 80 AU are visible during the initial and final epoch. The 
other runs show identical results. The lower panel of Fig. 6 shows in more detail the deviations of 
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Fig. 6. — : Accretion rate (upper panel) and deviations of the accretion rates from the simulation 
run with the smallest sink cell radius of r m ; n = 1 AU (lower panel) as a function of time for four 
different sizes of the spherical sink cell. 
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the simulations from the run with r m ; n = 1 AU. The simulation runs with r m j n — 10, 5, and 1 AU 
stay identical. In these simulations the dust sublimation radius, which can be roughly estimated 
to 20 to 30 AU for a corresponding 20 to 30 M star, is included in the computational domain 
before the onset of radiation pressure occurs at roughly 25 kyr. On the other hand, the largest sink 
cell of r min = 80 AU exceeds the dust sublimation radius r su bi- The resulting accretion rate of the 
corresponding run oscillates around the results from the more precise simulations with a maximum 
deviation of 10% mostly at the end of the simulation, when the radiation pressure starts to revert 
the accretion flow throughout the whole domain. Due to the fact that the deviations are oscillating 
and the fact that the strongest deviations occur at the end of the simulation where the accretion 
rate is already an order of magnitude lower than at the beginning, the four simulations yield the 
same final mass of the proto-star. Subsequent one-dimensional simulations presented make use of 
a radius of the central sink cell of r m ; n = 1 AU. 

4.2. Parameter scan of the initial pre-stellar core mass: The upper mass limit of 

spherically symmetric accretion 

4-2.1. Simulations 

The simulations presented so far were performed to fix the remaining free numerical parameters, 
namely, the grid resolution and the size of the central sink cell. We now study the collapse of 
massive pre-stellar cores for four different initial core masses ]\d COYC ranging from M corc — 60 M 
up to 480 Mq. The initial conditions and numerical parameters for these runs are described in 
Sect. 3. The simulations are performed with an inner boundary of the computational domain of 
Tmin = 1 AU and 128 grid cells with logarithmically increasing resolution towards the center. The 
size of the innermost grid cell of the computational domain is (Ar) m j n = 0.08 AU. We follow the 
evolution of the system until no mass is left in the computational domain. Part of this mass is 
accreted onto the central massive star and part is expelled over the outer boundary by radiative 
forces. The resulting accretion histories are displayed in Fig. 7 as a function of the actual stellar 
mass. 

4-2.2. Conclusions 

The mass and the luminosity of the forming massive star grow with time. The radiation 
pressure of the direct stellar irradiation as well as from the thermal infrared dust emission increase 
and become stronger than gravity ultimately. Therefore, the accretion rate drops down and the 
massive star has grown to its final mass. 

The individual force densities as as function of the radius through the spherically symmetric 
pre-stellar core are displayed at a snapshot in time, where the radiative force starts to trigger the 
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Table 4. Overview of initial conditions 



Symbol 


Value 


Quantity 


T 


20 K 


temperature of the pre-stellar core 


(M dust /M gas ) 


1% 


dust to gas mass ratio 


\u r \ 





radial velocity 


\ue\ 





polar velocity 


= \u<j>\/R 


5 * 10~ 13 Hz 


azimuthal angular velocity in 2D 


'"max 


0.1 pc 


outer radius of the pre-stellar core 


p{r) 


r -2 


density slope of the pre-stellar core 




60 to 480 M Q 


mass of the pre-stellar core 




10 20 30 40 

M* [M Q ] 

Fig. 7. — : Accretion rate M* as a function of the actual stellar mass M* for four different initial 
pre-stellar core masses of M core = 60 Mq up to 480 M . The spherically symmetric accretion 
models yield an upper mass limit of the final star of M^ D < 40 Mq. 
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stopping of the in-fall motion in Fig. 8. These forces are later on compared with the corresponding 
forces of the disk accretion models. 

The final star does not reach a mass higher than 40 M in any of the simulations. This 
limit is in good agreement with previous research studies. It lies in the allowed range of 25 - 
60 M Q determined by Larson & Starrfield (1971). Kahn (1974) predicted in his analytical study 
the formation of a 40 M Q star and Yorke & Kriigel (1977) formed a 36 M Q star in their radiation 
hydrodynamics simulation of a 150 M Q collapsing core. 

The oscillations of the accretion rate during the stopping of the in-fall motion are due to a 
negative feedback effect of the accretion luminosity: By increasing the initial mass of the pre- 
stellar core from 60 Mq up to 480 Mq, the amplitude of the accretion rate and therefore the 
accretion luminosity increases as well. Due to the resulting stronger radiative force, the increase of 
accretion luminosity leads to a de-acceleration of the accretion flow, which results in a reduction of 
the corresponding accretion luminosity. This negative feedback yields a highly episodic accretion 
history. The effect is stronger in cases, where the ratio of the accretion luminosity to the stellar 
luminosity is high, i.e. the effect is stronger for more massive cores. Such an oscillating phase was 
also previously detected in the simulations by Yorke & Kriigel (1977). 

The fact that the final mass of the star in the most massive case M core = 480 Mq is lower (M* ~ 
31 Mq) than for the cores that initially had less mass, should be taken with care: In simulations 
with such high oscillations, the influence of the underlying stellar evolution model increases strongly. 
To analyze the details of this time dependent interaction of the stellar evolution and the accretion 
flow, a self-consistent treatment of the proto-stellar's evolution and its environment should be 
considered. 

5. Disk accretion 

The most massive stars known cannot be formed by spherically symmetric accretion. As 
shown in the last section, the radiative forces in a spherically symmetric envelope lead to a cut-off 
of the accretion phase. The high luminosity of a massive star heats the region in its vicinity to 
such a high temperature that the resulting thermal radiation pressure overcomes the gravitational 
force. The radiation pressure stops, and finally reverses the accretion flow. Besides this theoretical 
issue, observations indicate the presence of angular momentum in all epochs of star formation, 
starting with the rotation of pre-stellar cores and finally resulting in rotating flattened circumstellar 
structures. Leaving perfectly spherical symmetry will thereby potentially help to overcome the 
radiation pressure problem. First, the presence of higher densities in the forming disk region 
results in a thinner shell, where the first absorption of stellar photons takes place. This enables 
an accretion flow to break through this region of direct stellar irradiation feedback more easily. 
Secondly, the feedback by radiation from dust grains, which actually stops the accretion in the 
spherically symmetric case, will be strongly reduced, because the majority of the radiative flux 
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Fig. 8. — : Snapshot of radial force densities (lower panel) and the density and temperature profile 
(upper panel) in the innermost core region taken from the collapse simulation of a M core = 120 M Q 
pre-stellar core at 20 kyr corresponding to a proto-stellar mass of about M* = 25 Mq. Due to the 
superior radiative force the spherically symmetric accretion models yield an upper mass limit of 
the final star of Mj D < 40 M . 
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from the irradiated inner rim of the disk will escape in the vertical direction through the optically 
thin disk atmosphere and therefore does not interact with the radially inward-streaming accretion 
flow. The different kinds of radiative feedback in spherical symmetry as well as in an axially 
symmetric disk geometry are illustrated in the final discussion Sect. 6. 

Analogously to the discussion of the spherical symmetric simulations, we present in the follow- 
ing the results of axially and midplane symmetric simulations of the collapse of rotating massive 
pre-stellar cores. Before being able to scan the parameter space of different initial core masses 
(Sect. 5.2), we determine the required resolution in convergence runs (Sect. A. 2) and fix the radius 
r min of the central sink cell in various simulations (Sect. 5.1). 

5.1. Parameter scan of the size of the sink cell: The influence of the dust 

sublimation front 

5.1.1. Simulations 

In the following, we study the influence of the radius r m ; n of the inner sink cell, which 
equals the inner computational boundary, on the resulting accretion rate onto the evolving mas- 
sive star. We check this dependency in four simulations with a radius of the inner sink cell of 
r min = 1) 5, 10, and 80 AU. The initial conditions and numerical parameters of these runs are 
described in Sect. 3 and the simulations are performed for an initial core mass of M core = 60 Mq. 
We follow the long-term evolution of the runs for at least 92 kyr. The resulting accretion rate onto 
the forming star as well as the mass growth of the central star are displayed in Fig. 9. 

5.1.2. Conclusions 

In the spherically symmetric models, we conclude that the numerical results do not depend 
on the radius r m [ n of the central sink cell unless it is smaller than the dust sublimation radius 
r subi from the point in time at which the radiative force overcomes gravity. These results cannot 
easily be transferred to the axially symmetric disk configuration. Including centrifugal forces, which 
compensate the gravity in the disk region, the chosen location r m i n of the inner boundary of the 
computational domain influences the resulting accretion rate in two distinguishable effects: 

Due to the fact that the circumstellar disk is growing in time from the inside outwards, a 
smaller sink cell leads to an earlier onset of the disk formation phase during the simulation. In 
other words, a fluid package with an initial position at (ri,9i) and an initial rotation of Oj hold 
a centrifugal radius of r ccnt = G j^^.^ sin (9i) with the included mass M(r«), see Eq. (12). If this 
centrifugal radius is smaller than the sink cell radius r m [ n , the fluid package is accreted onto the 
forming star during the so-called free fall epoch at the beginning of the simulation. This effect is 
associated with the gas physics (hydrodynamics) of the pre-stellar core, because the gas represents 
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t [kyr] 

Fig. 9. — : Stellar mass M* (upper panel) and accretion rate M* (lower panel) as a function of 
time t for different radii r m ; n of the central sink cell in the collapse simulation of a rotating 60 M Q 
pre-stellar core. 
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99% of the mass of the pre-stellar core. A second important effect depending on the chosen sink cell 
radius is related to the dust and therefore to the radiation physics. The region in the vicinity of the 
forming massive star will be heated up to temperatures beyond the dust sublimation temperature. 
Therefore a gap is formed between the central star and the dust disk. Under the assumption 
that the absorption by gas in this gap is much smaller than the absorption by dust grains behind 
the dust sublimation front the inner rim of the dust disk determines the region of the first stellar 
radiative impact onto the accretion flow. Also the most important radiative feedback due to thermal 
re-emission by dust grains sets in directly behind this irradiated heated region. Contrary to the 
spherically symmetric case, the high-density disk region bypasses most of the re-emitted radiation 
into the vertical direction, i.e. the total radiation field is composed of the isotropic stellar irradiation 
and a highly anisotropic (secondary) thermal radiation field. 

Fig. 9 illustrates clearly both, the mass and the radiative effect, related to an artificial inner 
cut-off of the gas and the dust disk respectively: As expected, the duration of the so-called free fall 
phase shortens with the radius r m ; n of the sink cell. This behaviour can fortunately be estimated 
analytically given the sink cell radius and the initial conditions of the pre-stellar core to account 
for the overestimation of the final mass of the forming star, if necessary. Moreover, this effect of 
the artificial inner rim of the gas disk results on the one hand in an overestimation of the final 
mass of the central star by approximately 1 M Q or below (upper panel of Fig. 9), but on the 
other hand influences the proceeding radiation hydrodynamic interactions in its environment only 
marginally (lower panel of Fig. 9). The corresponding accretion rates after the disk formation are 
not influenced at all. This result is quite reasonable keeping in mind that the balance of radiative 
and gravitational forces can be described in first order by the luminosity to mass ratio L to t/-^* 
of the central massive star, which only changes marginally with another choice of the size of the 
central sink cell. 

But the radiative impact due to the artificial cut-off of the inner dust disk regime for sink cell 
radii larger than the dust sublimation front has dramatic effects. In the case of r m ; n = 80 AU the 
artificial shift of the region of the dust radiation interactions terminates the short disk accretion 
phase abruptly, leading to a completely wrong evolution of the central star, the disk as well as the 
large scale environment. The reason for this dramatic change in the radiation physics is that the 
lower density region of the circumstellar disk at 80 AU is (in contrast to the real inner rim of the 
dust disk at roughly 20 AU) not opaque enough to generate a strong anisotropy of the thermal 
radiation field. Therefore the strong isotropic part of the thermal radiation field is able to stop the 
accretion analogous to the spherically symmetric flow calculations. 

Due to the importance of this inner core region for the associated interaction of the radiation 
with the accretion flow it seems to be unavoidable to include the whole dust disk down to its 
inner rim in the computational domain (cp. Fig. 10). This defines an upper limit of the radius 
r min of the central sink cell, which has to be smaller than the dust sublimation radius r su bi in 
the midplane from that point in time at which the radiative force has grown to a competetive 
magnitude compared to the viscous force driving the accretion flow. Subsequent simulations meet 
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this concern by using an adequate central sink cell radius of r m [ n = 10 AU. Otherwise, for an inner 
sink cell radius r m j n larger than the dust sublimation radius r su bi, the region of radiative feedback 
is artificially shifted to higher radii including a strong decrease in density, opacity, and gravity. The 
resulting strong heating of the disk region behind the radius r m j n > r su bi, which 'realistically' would 
be shielded from the stellar irradiation by the inner parts of the disk, leads to a diminishment of 
the shielding/shadowing property of the massive accretion disk. In case the radiation field retained 
therefore its isotropic character in major parts, the radiation pressure stops the emerging disk 
accretion phase, similar to the spherically symmetric case. 

This dependency of the radiation pressure on the radius of the sink cell explains also the 
abrupt end of the accretion phase in the simulations by Yorke & Sonnhalter (2002). They presented 
simulations of collapsing pre-stellar cores of M core = 30 M , 60 M Q and 120 M . The radius of 
their inner sink cell was chosen proportional to the initial mass of the core to be 40, 80, and 160 AU 
respectively. As shown in this parameter scan (cp. Fig. 9) such huge sink cells lead to an artificial 
cut-off of the dust disk and result therefore in unphysically strong radiative feedback. Therefore, 
we are definitely sure that this yields also the abrupt and early end of the accretion phase in the 
simulations by Yorke & Sonnhalter (2002). 

Contrary to our conclusion, Krumholz et al. (2009) stated that the much longer accretion 
phases in their own simulations compared to Yorke Sz Sonnhalter (2002) are a result of non-axially 
symmetric modes in the outflow region. The physical argument in that case is that their simulations 
show a so-called "3D radiative Rayleigh- Taylor instability" in the radiation pressure driven outflow, 
which results in further mass accretion onto the circumstellar disk from the bipolar direction instead 
of a steady outflow feature. In axially symmetric simulations the feeding of the disk would therefore 
not be possible. Our axially symmetric collapse simulations, presented in the following section, show 
a stable radiation pressure driven outflow and the forming circumstellar disk gains enough mass 
from the huge mass reservoir of the envelope to maintain its shielding property over several free 
fall times, in fact over a longer period ever simulated in previous research studies. 

5.2. Parameter scan of the initial pre-stellar core mass: Breaking through the 
upper mass limit of spherically symmetric accretion 

5.2.1. Simulations 

The spherically symmetric (one-dimensional) collapse simulations of massive pre-stellar cores 
yield a maximum stellar mass of less than 40 Mq independent of the initial core mass M core > 60 M Q 
due to radiative feedback. We attack this radiation pressure barrier in two-dimensional axially and 
midplane symmetric circumstellar disk geometry now. The implications of this change of geometries 
are illustrated at full length in the final discussion Sect. 6. We performed four simulations with 
the default initial conditions described in Sect. 3 and the fixed numerical parameters presented in 
Sect. A. 2 and 5.1, namely a maximum resolution of 1.27 AU x 1.04 AU, and an inner sink cell radius 




Fig. 10. — : Resolving the dust sublimation front. The image shows a snapshot at 100 kyr after 
start of the collapse of a 120 M pre-stellar core. 

Color: Gas density from 10 -20 up to 2 * 10~ 12 g cm -3 in logarithmic scale. 
Contour-lines: Dust to gas mass ratio from up to 1% in linear scale. 

The zoom-in illustrates the fact that the dust sublimation front is resolved smoothly over several 
(here roughly 8) grid cells in the radial direction. 
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of 10 AU. The different initial core masses of M core = 60 M , 120 M , 240 M , and 480 M are 
chosen analog to the scan of the initial core mass parameter in the spherically symmetric case. The 
resulting accretion histories as a function of the actual stellar mass are displayed in Fig. 11. 
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Fig. 11. — : Accretion rate M* as a function of the actual stellar mass M* for four different initial 
core masses M core = 60 M , 120 M , 240 M , and 480 M . Also the periods of accretion are 
mentioned for each run. The vertical line marks the upper mass limit found in the spherically 
symmetric accretion models. The collapse models of slowly rotating pre-stellar cores clearly break 
through this upper mass limit of the final star of M* D < 40 M . 



5.2.2. Conclusions 

As expected, the lowest mass case of M core = 60 M results finally in a less massive central 
star than the corresponding run in spherical symmetry (without rotation) simply due to the fact 
that the additional angular momentum results in centrifugal forces, which counteracts the accretion 
flow driven by gravity and viscosity. In face of this additional centrifugal forces, for higher mass 
pre-stellar cores the slowed down accretion flux breaks easily through the upper mass limit of the 
final star of Af^ D < 40 M found in the spherically symmetric accretion models! 

The reason for that breakthrough can be displayed by a closer look at the driving force densities 
in the evolved pre-stellar core, plotted in Figs. 12 to 17. All figures represent a snapshot of the 
■Mcore — 120 M case at 60 kyr after start of the simulation. At this point in time, the actual mass 
of the central massive star is roughly 40 M , representing the spherically symmetric upper mass 
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limit found in previous simulations (Sect. 4.2). In contrast to the spherically symmetric models, 
the geometry of the proto-stellar environment can now be divided into a very dense circumstellar 
disk and the lower density envelope. We visualized exemplary the actual density, velocity, and the 
acting forces in the radial direction for both regimes, Figs. 12 to 14 for the midplane of the accretion 
disk, Figs. 15 to 17 for a polar angle of 30° above the midplane. In the midplane the gravity and 
centrifugal force are one to two orders of magnitude higher than the thermal pressure and up to 
three orders of magnitude higher than the radiative and viscous force. The upper panel of Fig. 12 
shows three individual regions of the midplane layer, in between the sign of the total force density 
changes. The gravity dominates the individual forces for the outer core regions (above 3000 AU) 
leading to a steady accretion flow onto the inner core region (Figs. 12 and 13). In the very inner 
part of the core around the massive star (below 200 AU) the gravity is balanced by the centrifugal 
force and in small part by the thermal pressure (Fig. 13). In this region, which we will refer to as 
the disk region hereafter, the shear viscosity yields a quasi-stationary accretion flow through the 
disk, which clearly exceeds the radiative force (Fig. 14). In between this disk region (< 200 AU) 
and the global in- fall region (> 3000 AU) the mass flux describes transient oscillations, because 
gravity, centrifugal forces and thermal pressure are not in equilibrium yet, as it is the case for the 
mass finally arriving the disk region. Although the total force density, displayed in the upper panel 
of Fig. 12 is directed in the outward direction between 200 and 3000 AU, the mass in this region is 
still in an inward motion (cp. Fig. 12, lower panel), i.e. the mass flow through the pre-stellar core 
is feeding the circumstellar accretion disk. The viscous force in the accretion disk is able to drive 
a steady accretion flow towards the evolving massive star of 40 M , because the radiative force 
is one to two orders of magnitude lower in this dense disk region than in the low density envelope 
(cp. Figs. 14 and 17). Observations of such a large scale flattened structure with a potentially 
embedded small scale accretion disk are i.e. described in Fallscheer et al. (2009) and Beuther et al. 
(2009). 

At an polar angle of 30° above the midplane this strong radiative force already accelerates the 
remnant mass in the radially outward direction through mostly the entire pre-stellar core (Fig. 16). 
Only at the outer rim of the core the previous in-fall motion is still visible. This distribution of the 
individual force densities confirms in high detail the expected result of the anisotropy of the thermal 
radiation field. Most of the radiative flux from the irradiated inner rim of the disk is bypassed in the 
vertical direction through the optically thin atmosphere of the circumstellar disk. Meanwhile, the 
accretion flow is reduced compared to the one-dimensional gravitational in-fall to a steady stream 
driven by the viscous properties of the accretion disk. In the envelope region of the pre-stellar 
core the radiative force reverts the in-fall motion and depletes the stellar surrounding similar to 
the spherically symmetric accretion models (cp. the corresponding individual force distributions in 
Figs. 15 and 8). 

In these axially and midplane symmetric disk accretion models no upper mass limit of the 
final star is detected so far, but the star formation efficiency declines for higher mass cores. The 
proceeding depletion of the envelope by radiative forces finally leads to a decrease of the density 
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in the midplane and the disk looses its shielding property. Without this shielding the radiation 
pressure starts to accelerate the remnant material in the outward direction. 



6. Discussion 

6.1. Radiation pressure feedback in a nutshell 

Since the isotropic and anisotropic feature of the spherical and disk accretion scenario at 
the dust sublimation front is of such a great importance in our simulation results, we illustrate 
these key attributes in more detail: In this research study, we investigated the influence of the 
stellar environment onto the radiation pressure problem in the formation of massive stars. We 
studied the accretion flow onto a high-mass star in a monolithic pre-stellar core collapse picture, 
as recommended by Whitney (2005) and McKee & Ostriker (2007). Under this assumption, the 
theoretical description of the accretion process onto a massive star has to deal with the interaction 
between the exerted radiation by the forming star with the accretion flow of gas and dust (Shu ct al. 
1987). In a perfectly spherically symmetric collapse, this interaction potentially stops the accretion 
onto the star entirely. In the static limit, radiation pressure overcomes gravity at the so-called 
generalized Eddington barrier 

where L*, M* and ft* denote the stellar luminosity, the stellar mass, and the dust opacity respec- 
tively, G is the gravitational constant and c is the speed of light. But the collapse of a pre-stellar 
core is far from being a static problem. The momentum transfer from the absorbed photons first 
has to slow down the in-falling envelope. For simplification purposes, we can divide the radiation 
pressure feedback into two components, as illustrated in Fig. 18. The first exchange of momentum 
takes place when the irradiation from the massive star is absorbed by the dust grains of the sur- 
rounding, i.e. behind the dust sublimation radius. The strongest force will thereby be produced 
by photons with shorter wavelengths, because they have a higher absorption probability and are 
more energetic. We will call this first interaction "UV feedback" abbreviated, although the fre- 
quency dependence of the broad stellar black body spectrum is clearly not negligible. Afterwards, 
these heated regions emit most of the photons at the dust temperature, yielding a much longer 
wavelength and a much longer mean free path than the direct stellar light. The interaction of this 
radiation with the enclosed gas and dust is therefore referred to as "IR feedback" . Our spherically 
symmetric collapse simulations confirm the outcome of previous studies that it is essentially the 
IR feedback that stops the accretion flow onto the forming star in spherical symmetry. Although 
each IR photon transfers less momentum to the dust than the highly energetic stellar UV pho- 
tons, the thermal dust emission acts onto the accretion flow on a much larger volume than the 
spatially confined absorption region of the stellar irradiation. Additionally the optical depth of 
the envelope decreases towards longer wavelength, so the IR feedback conteracts the accretion flow 
in the outer core regions yielding less gravity. Different approaches to overcome this barrier for 
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Fig. 12. — : Total force density |/tot( r )| (upper panel) as well as density p(r) and radial velocity 
u r (r) (lower panel) as a function of radius r through the disk's midplane. The snapshot was taken 
at 60 kyr after start of the simulation, corresponding to a central stellar mass of roughly 40 M . 
The individual force densities along this line of sight through the total and the inner core region 
are displayed in Figs. 13 and 14. 
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Fig. 13. — : Gravity, centrifugal, and thermal pressure force as a function of radius through the 
disk's midplane. The snapshot was taken at 60 kyr after start of the simulation, corresponding to 
a central stellar mass of roughly 40 Mq. The radiative and viscous forces are orders of magnitude 
smaller than the illustrated ones, but become important in the inner disk region, where the stronger 
forces are in equilibrium. The radiative and viscous force densities along this line of sight through 
the inner core region are displayed in Fig. 14. 
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Fig. 14. — : Viscous and radiative force density of the inner core region as a function of radius 
through the disk's midplane. The snapshot was taken at 60 kyr after start of the simulation, 
corresponding to a central stellar mass of roughly 40 Mq. 
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Fig. 15. — : Total force density | /tot (01 (upper panel) as well as density p(r) and radial velocity 
u r (r) (lower panel) as a function of radius r at 30° above the disk's midplane. The snapshot was 
taken at 60 kyr after start of the simulation, corresponding to a central stellar mass of roughly 
40 M.q. The individual force densities along this line of sight through the total and the inner core 
region are displayed in Figs. 16 and 17. 
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Fig. 16. — : Gravity, centrifugal, thermal pressure, and radiative forces as a function of radius at 
30° above the disk's midplane. The snapshot was taken at 60 kyr after start of the simulation, 
corresponding to a central stellar mass of roughly 40 Mq. The individual force densities along this 
line of sight through the inner core region are displayed in Fig. 17. 
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Fig. 17. — : Gravity, centrifugal, thermal pressure, radiative, and viscous force density of the inner 
core region as a function of radius at 30° above the disk's midplane. The snapshot was taken at 
60 kyr after start of the simulation, corresponding to a central stellar mass of roughly 40 M . 
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Fig. 18. — : Schematic view of the radiative forces onto the accretion flow in spherical symmetry. 
The radiative feedback is divided into direct stellar irradiation and secondary re-emitted photons 
by dust grains. 



spherically symmetric accretion flows onto massive stars were considered in the past. The general- 
ized Eddington barrier depends only on the stellar evolution (L*/M*) and on the dust properties 
(k*). Wolfire & Cassinelli (1987) studied the necessary change of dust properties to enable further 
accretion, but the restrictions they derived seem to be unrealistic. 

Without a doubt, star formation is rarely a perfectly spherically symmetric problem. Initial 
angular momentum of the collapsing pre-stellar core leads to the formation of a circumstellar disk as 
well as polar cavities. Compared to the case of spherically symmetric accretion, the disk geometry 
changes the radiation pressure feedback dramatically, see Fig. 19. Going from a spherically sym- 
metric in-fall to an axially symmetric disk geometry can help to overcome both - the UV as well as 
the IR - radiation pressure feedback: Developing radiation hydrodynamical (Klahr & Bodenheimer 
2003), magneto-rotational (Balbus & Hawley 1991; Hawley & Balbus 1991; Balbus 2003) and self- 
gravitating instabilities (Yang et al. 1991; Laughlin & Bodenheimer 1994; Bodenheimer 1995) in 
the accretion disk will transfer angular momentum to outer radii allowing a steady mass accretion 
radially inward. The additional ram pressure in the radiatively shielded parts of the disk will pos- 
sibly push the mass over the thin shell of the UV feedback (Nakano 1989). Secondly, and most 
important, the irradiated and therefore heated regions of the disk will mainly cool in the vertical 
direction through the optically thin disk's atmosphere, strongly restraining the IR radiation pres- 
sure in the radial direction. If the latter process occurs at the innermost part of the disk so that the 
radiation can escape directly through the bipolar cavity, this effect is also known as the so-called 
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Fig. 19. — : Schematic view of the "UV"- and "IR" -component of the radiation pressure acting in 
an axially symmetric circumstellar disk geometry. 

'flashlight-effect' (Yorke & Sonnhalter 2002; Krumholz et al. 2005). 

The interaction of the radiation with the accretion flow is very sensitive to the numerical 
treatment of radiation transport. The FLD approximation, which is a standard technique in mod- 
ern radiation hydrodynamics codes for astrophysical fluid flows, fails to compute the correct flux 
between the first transition region from the dust depleted zone around the massive star and the op- 
tically thick disk leading to an incorrect temperature distribution in the irradiated regions (see e.g. 
Yorke & Kriigel 1977; Boley et al. 2007). Also simplifying the stellar black body spectrum by using 
frequency averaged Planck mean opacities leads to a thinner shell of the direct stellar irradiation 
feedback and a stronger heating of the corresponding dust, which afterwards yields a higher IR 
feedback. Hence, accounting for the frequency dependence of the stellar spectrum seems to be a 
crucial point. 

The most violent interaction of the stellar irradiation with the accretion flow takes place at 
and directly behind the first absorption peak. The location of the first absorption layer is repre- 
sented by the dust sublimation front, where the local dust temperature falls below the evaporation 
temperature of the dust grains. A systematic study of the radiation pressure feedback on the for- 
mation of massive stars therefore implies the need to resolve the ongoing radiation and accretion 
physics down to the dust sublimation front. A formation of massive stars by breaking through the 
ionization boundary into regions of sublimated dust grains was studied for spherically symmetric 
accretion flows (Keto 2003) as well as for two-dimensional effects in the so-called small radius limit 
(Jijina & Adams 1996). Aside from the important contribution of the proceeding physics at the 
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dust sublimation front, no previous non-spherically symmetric numerical research has been done 
so far, presumably due to resolution issues. 

6.2. Comparison to previous numerical research in the field 

In our simulations the star in the center of the accretion disk grows far beyond the upper 
mass limit found in the case of spherical accretion. Indeed, the final massive stars are the most 
massive stars ever formed in a multi-dimensional radiation hydrodynamics simulation so far. The 
quantitative final results as well as a comparison to previous multi-dimensional radiation hydrody- 
namics studies is presented in Table 5. The research studies clearly differ in the evolutionary 
time simulated. We improved this by roughly an order of magnitude. The 20 kyr and 75 kyr of 
evolution in Krumholz et al. (2007) and Krumholz et al. (2009) represent approximately one-third 
and slightly more than one free fall time of the pre-stellar core respectively. Despite the frequency 
dependent radiation transport and the high resolution down to 1.27 AU in our simulations, we are 
able to follow the evolution of the accreting system up to 14, 10, 7, and 2 free fall times for an initial 
core mass of 60, 120, 240, and 480 M Q respectively, including the whole stellar accretion phase. To 
state this clearly, this is only possible due to our self-restriction to axial symmetry in these runs. 
In the simulations by M. Krumholz further accretion seems to be the natural continuation of the 
runs. Our simulation series of the disk accretion scenario with varying initial core masses shows 
a decrease of the star formation efficiency towards higher mass cores as a result of the growing 
radiation pressure feedback. 

Simultaneously to the bypass of the thermal radiation by the massive accretion disk, a stable 
wide-angle bipolar outflow with velocities of the order of 100 km s" 1 is launched by the radiation 
pressure. In these axially symmetric simulations, we did not detect any evidence for a radiative 
instability of these outflow regions, like observed in the simulations by Krumholz et al. (2009). 
Following the explanatory notes in Krumholz et al. (2009), this is due to the fact, that this insta- 
bility requires non-axial symmetric modes to occur. A detailed study of this regime seems to be 
necessary to clearly understand the underlying physics of this requirement. On the other hand, 
our simulations show a strong release of radiation pressure in the bipolar direction, growing in 
angle with time, which fits to the observed broadening of outflows in massive star forming regions 
(Beuther & Shepherd 2005). But to state this clearly, we are sure that a complete description of the 
jet or outflow physics cannot be done without taking care of the dominant magnetic field effects. 

Fig. 20 shows the time dependent fractions of masses, divided into the mass inside the com- 
putational domain, the mass of the forming star (from the flux over the inner boundary), and the 
mass loss by the radiation pressure driven outflow (from the flux over the outer boundary). 

The different epochs of the collapse of the 120 M pre-stellar core are illustrated in Fig. 21. 
The subfigures display the initial condition, the disk formation and evolution, the outflow launching, 
and the end of the accretion phase in several snapshots of the density structure. The corresponding 
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Table 5. Overview of multi-dimensional radiation hydrodynamics simulations of massive star 

formation 





-'"core 


^end 


M* 


SFE 


Authors 


(M ) 


(kyr) 


(M ) 


(%) 




30 


25 


31.6 




Yorke & Sonnhalter (2002) 


60 


45 


33.6 






120 


70 


42.9 






100 (A) 


20+ 


5.4 (+ 3.4) 


> 5.4 


Krumholz et al. (2007) 


100 (B) 


20+ 


8.9 (+ 2.4) 


> 8.9 




200 


20+ 


8.6 (+ 6) 


> 8.6 


I. (2( )9) 


100 


75+ 


41.5 + 29.2 (+ 28.3) 


> 70.7 




60 


939 


28.2 


47.0 




120 


489 


56.5 


47.1 


This study 


240 




92.6 




226 


38.5 




480 


41+ 


137.2 (+ 67.8) 


> 28.5 but < 42.7 



Note. The columns from left to right state the authors, the initial pre-stellar core mass, 
the evolutionary time simulated, the final star mass, as well as the corresponding star formation 
efficiency. A "+" in the t e nd column means that the accretion phase is not simulated until the end 
yet. In the case of the simulation by M. Krumholz, only the formation of the most massive stars 
are considered here; all other stars formed have masses below 1 M . In the case of Krumholz et al. 
(2007) the "(A)" and "(B)" in the M core column mark the usage of different perturbation fields of 
the initial state (same labels as in the original paper) and the M* column gives additionally the 
remnant disk mass around the primary star. In the case of Krumholz et al. (2009) and the 480 M 
case of our own study the M* column gives additionally the remnant disk plus envelope mass. In 
the case of Yorke & Sonnhalter (2002) only the simulations with frequency dependent radiation 
transport are considered. 
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Fig. 20. — : Evolution of mass flows over the inner and outer computational boundary. 

Solid line: Mass inside the computational domain. 

Dashed line: Mass of the forming massive star in the center. 

Dot-dashed line: Mass loss due to the radiation pressure in the radially outward direction. 
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animation is available in the online-materials. 

6.3. Limitations of our approach 

A minor flaw of our studies of the radiation hydrodynamics around the dust sublimation front 
is the usage of a simple constant gas opacity of Atg as — 0.01 cm g for the completely evaporated 
regions around the forming star. In other words, the neighborhood of the star remains optically 
thin for the stellar irradiation up to the dust sublimation front. From our results of the parameter 
scan of the inner sink cell radius, we conclude that the prior absorption of the stellar irradiation 
in dust-free, but potentially optically thick regions, would even enhance the crucial anisotropy of 
the radiation field detected in our simulations. In this sense, the ignorance of the detailed optical 
properties of the gas phase does not imply any loss of generality. 

The usage of a central sink cell of a specific radius implies further assumptions, e.g. we assume 
that the mass flow into the central sink cell is accreted by the central star. If, e.g. the mass flow 
is transferred into an outflow or jet in this inner region near the stellar surface, the stellar growth 
would be decreased and the final masses of the stars given here represent upper mass limits. 

Due to the fact that higher accretion rates lead first to a delay of the star's approach to the 
main sequence and secondly to a higher ratio of accretion luminosity to stellar luminosity (which 
increases the importance of the correct knowledge of the stellar radius), the details of the treatment 
of stellar evolution become important especially in the case of the highest mass cores studied. An 
even more realistic approach than the usage of tabulated tracks could be achieved by including a 
stellar evolution code such as Hosokawa & Omukai (2009) to calculate the ongoing stellar physics 
in the sink cell consistently in time. 

To mimic the effect of angular momentum transport by evolving instabilities in the accretion 
disk, we made use of the a-viscosity model by Shakura & Sunyaev (1973). Nevertheless, all our code 
development (radiation transport and self-gravity) and model setup is already in three-dimensional 
formulation. Hence we started to expand our simulations into full 3D. First, this will allow to com- 
pute the angular momentum transport (by gravitational torques) consistently with the formation 
and evolution of the accretion disk. Secondly, although we will not study the fragmentation of the 
outer core regions, the stability or potential fragmentation of the forming massive circumstellar 
disk can be addressed. 

7. Summary 

We performed high-resolution radiation hydrodynamics simulations of monolithic pre-stellar 
core collapses including frequency dependent radiative feedback. A broad parameter space of 
various numerical configurations and initial conditions was scanned. The dust sublimation front in 
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Fig. 21. — : Simulation snapshots from a collapse of a 120 M pre-stellar core. All images show the 
same scale of the whole pre-stellar core with an initial diameter of 0.2 pc. 
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Fig. 21. — : Continuation of Fig. 21 
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the vicinity of the forming star could be resolved down to 1.27 AU. The evolution of the system 
was computed over its whole accretion phase of several 10 5 yr. The usage of frequency dependent 
ray-tracing in our newly developed radiation module denotes the most realistic radiation transport 
method used in multi-dimensional hydrodynamic simulations of massive star formation by now. 
The broad parameter studies, especially regarding the size of the sink cell and the initial core mass, 
reveal new insights of the radiative feedback onto the accretion flow during the formation of a 
massive star: 

In the case of spherically symmetric accretion flows, we confirm the results of previous research 
studies (Larson & Starrfield 1971; Kahn 1974; Yorke & Kriigel 1977; Wolfire & Cassinelli 1987) that 
the thermal radiation pressure by re-emitted photons behind the dust sublimation front overcomes 
gravity, stops the accretion flow, and finally reverts the in-falling envelope. The upper mass limit 
of spherically symmetric accretion for our specific dust (Laor & Draine 1993) and stellar evolution 
model (Hosokawa & Omukai 2009) constrains the final stellar mass to be less than 40 M Q . 

In the case of disk accretion, the thermal radiation field generates a strong anisotropic feature, 
similar to the flashlight effect, which focus lies on the escape of radiation through optically thin 
cavities. We found that it is strict necessary to include the dust sublimation front in the computa- 
tional domain, to reveal the persistent anisotropy during the long-term evolution of the accretion 
disk. This requirement as well as a steady feeding of the accretion disk from the outer core regions 
maintain the anisotropic structure of the thermal radiation field. The short accretion phases of 
the disks in the simulations by Yorke & Sonnhalter (2002) are a result of the fact that they did 
not include the dust sublimation front in their simulations, as clearly shown in our result of the 
parameter scan of the size of the central sink cell (see Sect. 5.1). Additional feeding of the disk 
by unstable outflow regions, as stated in Krumholz et al. (2009), would enhance this anisotropy 
but is not necessary. As a consequence, we conclude that the radiation pressure problem in the 
formation of massive stars can be reduced to the question if the non-spherically symmetric stellar 
environment is dense or opaque enough to generate a strong anisotropy of the thermal radiation 
field. 

These mechanisms allow the central star to increase its mass far beyond the upper mass 
limit found in the case of spherical accretion flows. For an initial mass of the pre-stellar host 
core of 60, 120, 240, and 480 M the masses of the final stars formed in our simulations of the 
disk accretion scenario add up to 28.2, 56.5, 92.6, and at least 137.2 M respectively. Indeed, 
the final massive stars are the most massive stars ever formed in a multi-dimensional radiation 
hydrodynamics simulation so far. 

This research has been supported by the International Max-Planck Research School for Astron- 
omy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). Author H. Klahr has been 
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gruppe 759 "The Formation of Planets: The Critical First Growth Phase" . Author H. Beuther ac- 
knowledges financial support by the Emmy-Noether-Programm of the DFG through grant BE2578. 
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A. Parameter scans of the resolution 

Numerical hydrodynamics simulations involve a discretization of the underlying equations of 
hydrodynamics given in continuous space, cp. Eqs. (2)-(4). This causes a discretization error, 
which in general vanishes for infinitely high resolution of the numerical solver method. To compute 
a specific quantity, such as the accretion history, with a specific accuracy therefore needs a specific 
resolution, which is necessary to damp the discretization errors down to the requested accuracy. In 
this way it is possible to guarantee the achievement of a converged result. Although this procedure 
is a must in numerical research to achieve reliable results, the overhead of cpu time needed for 
convergence runs inhibits their realization in most present-day astronomical simulations, especially 
in multi-dimensional radiation hydrodynamics, which are performed almost at the upper limit of the 
computational power of the available clusters. To fix the number of grid cells, which are necessary 
for a correct representation of the radiation fluid interactions, we perform several simulations with 
varying resolution. Focusing on the inner regions of the pre-stellar core, the radial cell sizes of the 
grid thereby grow logarithmically from inside out as described in Sect. 2.1. 

A.l. ID convergence runs 

A. 1.1. Simulati ons 

The initial conditions and numerical parameters of these convergence runs are described in 
Sect. 3. The simulations are performed for an initial core mass of M COTC = 60 M and with an inner 
sink cell radius of r m j n = 1 AU. We follow the long-term evolution of the system for at least 163 
kyr, representing 2.4 free fall times of the pre-stellar core. The resulting mass growth M*(t) of the 
centrally forming star is displayed in Fig. 22. 

A. 1.2. Conclusions 

The lowest resolution run {N r = 32) fails to compute the correct amount of accretion already 
during the mostly isothermal initial free fall phase (up to 25 kyr). For higher resolution runs 
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Fig. 22. — : Stellar mass M* as a function of time t for four different resolutions of the spherically 
symmetric pre-stellar core collapse simulations. The number of grid cells N r varies from 32 to 256, 
corresponding to a size of the smallest grid cell of (Ar) m i n = 0.36 AU to 0.04 AU respectively. 
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with 64 grid cells or more, the mass growth of the forming star is identical during this phase. At 
a later evolutionary epoch, when radiative feedback becomes important, simulations with higher 
resolution lead generally to a slower mass growth. The deviation of a specific run to the next one 
with double resolution shrinks towards higher resolution, so the simulations fulfill the requirement 
of a monotonic convergence. Our one-dimensional simulations with varying initial core masses 
(Sect. 4.2) use 128 grid cells in the radial direction and an inner radial boundary at r m ; n = 1 AU 
corresponding to a grid size of (Ar) m ; n = 0.08 AU for the innermost grid cell. 

A. 2. 2D convergence runs 

A. 2.1. Simulations 

To fix the number of grid cells necessary for computing the correct physics of the radiation 
fluid interaction, we performed several simulations with varying resolution in the two-dimensional 
setup, too. The basic initial conditions and numerical parameter used for these convergence runs 
are described in Sect. 3. The simulations are performed for a core mass of M core = 60 M and the 
inner boundary of the computational domain is located at r m ; n = 10 AU. We followed the evolution 
of the collapsing core up to 33 kyr (0.5 free fall times) for the highest resolution case yet and up 
to several hundred kyr (about 10 free fall times) for a long-term convergence run. The accretion 
history and the corresponding mass growth of the centrally forming star are displayed in Figs. 23 
and 24. 

A. 2. 2. Conclusions 

In contrast to the purely one-dimensional in-fall (Sect. A.l), the centrifugal forces slow down 
the radially proceeding dynamics. So the usage of 64 grid cells in the radial direction, corresponding 
to a radial grid size of the innermost cells of (Ar) m ; n = 1.27 AU, yields a converged result for the 
accretion rate onto the forming high-mass star. The low-resolution run with only 32 grid cells in the 
radial direction clearly fails to compute the correct onset of disk formation after 8 kyr. Due to the 
clear dominance of the motion of gas in the radial direction during the initial 'free fall' phase up to 
roughly 8 kyr the accretion rates of this epoch are independent of the resolution used in the polar 
direction. The required resolution in the polar direction to compute a converged result also during 
later epochs remains notably poor, reflecting the fact that the complex radiation hydrodynamics 
interactions act mostly in the radial direction. This result confirms the expedient choice of spherical 
coordinates in monolithic core collapse simulations. Higher resolution of the polar stratification of 
the forming circumstellar disk mostly influences the cooling of the irradiated and viscously heated 
midplane layer. The usage of only 4 or 8 grid cells in the polar direction therefore results in stronger 
fluctuations of the accretion flow, which vanish in the higher resolution runs (clearly visible in the 
lower panel of Fig. 23). On the other hand the runs with low resolution in the polar direction 
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Fig. 23. — : Stellar mass M* (upper panel) and accretion rate M* (lower panel) as a function of 
time t for five different resolution to determine the adequate number of grid cells necessary for the 
collapse simulations of the rotating pre-stellar cores. 
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Fig. 24. — : Stellar mass M* (upper panel) and accretion rate M* (lower panel) as a function of time 
t for two different resolution in a long-term convergence study up to the end of the disk accretion 
phase. 
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underestimate the mass growth of the forming star only slightly (upper panel in Fig. 23). The 
deviations of each run to the next run in higher resolution shrink towards higher resolution, that 
means the simulation series yields a monotonous convergence. The long-term convergence study 
(Fig. 24) clearly shows that the point in time when the disk looses its shadowing property depends 
on the polar resolution of the circumstellar disk. Higher resolution of the disk's stratification results 
in a stronger anisotropy of the thermal radiation field and therefore minimizes the radiation pressure 
on the accretion flow. 

The runs with 64 x 16 and 128 x 32 grid cells show fully converged results even during the epoch 
of the most rapid changes at the onset of disk formation at 8 kyr. The spike in the accretion rate 
downwards during this onset represents the short period in time, in which for the first time a fluid 
package from the outer core region arrived at the innermost radius r m ; n with high enough angular 
momentum to compensate the stellar gravity. Quickly hereafter the following mass builds up a 
circumstellar disk, in which the shear viscosity yields an angular momentum transfer outwards 
resulting in a steady accretion rate anew. At later evolutionary phases the amplitude of the 
accretion rate is mostly a result of a quasi-stationary accretion flow inwards and an interactive 
radiative flux in the outward direction, which smoothly grows proportional to the luminosity of 
the forming massive star. The deviations of the individual runs during this more evolved and 'less 
violent' phase shrink again for all resolutions studied. Our two-dimensional simulations presented 
use 64 x 16 grid cells as the default setup corresponding to a grid size of (Ar x rA9) m i n = 1.27 AU 
x 1.04 AU for the innermost grid cells. 
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